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ABSTRACT 

We present numerical results on three-dimensional (3D) hydrodynamic core-collapse simulations of an 
11.2M© star. By comparing one-(lD) and two-dimensional(2D) results with those of 3D, we study how the 
increasing spacial multi-dimensionality affects the postbounce supernova dynamics. The calculations were per- 
formed with an energy-dependent treatment of the neutrino transport that is solved by the isotropic diffusion 
source approximation scheme. In agreement with previous study, our ID model does not produce explosions 
for the 11.2 M Q star, while the neutrino-driven revival of the stalled bounce shock is obtained both in the 2D 
and 3D models. The standing accretion-shock instability (SASI) is observed in the 3D models, in which the 
dominant mode of the SASI is bipolar {t = 2) with its saturation amplitudes being slightly smaller than 2D. 
By performing a tracer-particle analysis, we show that the maximum residency time of material in the gain 
region becomes longer in 3D due to non-axisymmetric flow motions than in 2D, which is one of advantageous 
aspects of 3D models to obtain neutrino-driven explosions. Our results show that convective matter motions 
below the gain radius become much more violent in 3D than in 2D, making the neutrino luminosity larger 
for 3D. Nevertheless the emitted neutrino energies are made smaller due to the enhanced cooling. Our results 
indicate whether these advantages for driving 3D explosions could or could not overwhelm the disadvantages 
is sensitive to the employed numerical resolutions. An encouraging finding is that the shock expansion tends 
to become more energetic for models with finer resolutions. To draw a robust conclusion, 3D simulations 
with much more higher numerical resolutions and also with more advanced treatment of neutrino transport as 
well as of gravity are needed, which could be hopefully practicable by utilizing forthcoming Petaflops-class 
supercomputers. 

Subject headings: supernovae: collapse — neutrinos — hydrodynamics 



1. INTRODUCTION 

Core-collapse supernovae have long drawn the attention of 
astrophysicists because they have many aspects playing im- 
portant roles in astrophysics. They are the mother of neu- 
tron stars and black holes; they play an important role for ac- 
celeration of cosmic rays; they influence galactic dynamics 
triggering further star formation; they are gigantic emitters 
of neutrinos and gravitational waves. They are also a major 
site for nucleosynthesis, so, naturally, any attempt to address 
human origins may need to begin with an understanding of 
core-collapse supernovae. 

Eve r since the first numerical simulation dColgate & White! 
119661) . the neutrino-heating mechanism, in which a stalled 
bounce shock is rev ived by neut r ino energy deposition 
to trigger ex plosions dWilsonl 119851: iBethe & Wilsonl Il985t 
iBethd Tl990l) . has been the working hypothesis of su- 
pernova theorists for these ~ 45 years. However, one 
impor ta nt lesson we hav e learn e d from | Rampp & Jankal 
d2000b: iLiebendorferet al] d200ll) : iThompson et all d2003l) ; 
ISumivoshi et al.l d2005l) who implemented the best input 
physics and numerics to date, is that the mechanism fails to 
blow up canonical massive stars in spherical symmetric (ID) 
simulations. Pushed by mo unting supernova o b servations of 
the bla s t morphology (e.g.. [Wang et al.l d2001l) ; iMaeda et al.l 
( 120081) : iTanaka et all d2009l) . and references therein), it is 
now almost certain that the breaking of the spherical sym- 
metry is the key to solve the supernova problem. So 
far a number of multidimensional (multi-D) hydrodynamic 



simulations have shown that hydrody namic motions asso 



Herant et all (fl994l); 
d!996l) : iFrveretalJ 



dated with convectiv e overturn (e.g. 
iBurrows etall d 19951) : Uanka & Miillei 
d2002l) ; iFrver] j2004 ) and the Stand i ng-Ac c retion-Shock 
Instab il ity (S A SI. e.g.. iBlond netai] (12003b: fSchecketal 
([2004] 12005): lOhnishi et al.l (120061 120071): iFoglizzo et al 
d2006l) : Ifwakami et all (I2008L 120091): lMurphv& Burrows 
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d2008l) : iFernandez & Thompsonl d2009bllal) . and references 
therein) can help the onset of the neutrino-driven explosion. 

In fact, the neutrino-driven explosions have been obtained 
in the following state-of -the-art t wo-dim ensional (2D) sim- 
ulations (e.g., table 1 in iKotakd d2011h ). Using the MuD- 
BaTH code which includes o ne of the best ava ilable neu- 
trino transfer approximations, iBuras et ail d2006l) firstly re- 
ported explos ions for a no n -rotati ng low-mass (11.2M©) 
progenitor of IWooslev et al.l d2002l). an d then for a 15M© 
progenitor of IWooslev & Weaver! ( 19951) with a moderately 
rapid rotation imposed dMarek & Jankal 120091) . By imple- 
menting a multi-group fl ux-limited diffusio n algorithm to 
the CH IMERA code (e.g. jBruenn et alJl2010h , lYakunin et ail 
d2010l) obtained expl osions for a non-rotating progenitors of 
IWooslev etail d2002l) in th e mass range of 12M and 25M Q . 
More recently, ISuwa et afl (120101) pointed out that a stronger 
explo sion is obtained for a r a pidly rotating 13M Q progeni- 
tor of iNomoto & Mashimotol d 19881) compared to the corre- 
sponding non-rotating model, in which the isotropic diffusion 
sourc e approximation (IPSA) for the spectral neutrino trans- 
port dLiebendorfer et af] 12009b is implemented in the ZEUS 
code. 
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This success, however, is opening further new questions. 
First of all, the explosion energies obtained in these 2D simu- 
lations are typically underpowered by one or two orders of 
magnitudes to explain the canonical supernova kinetic en- 
ergy (~ 10 51 erg). Moreov er, the softer nuclear equat ion of 
state (EOS), such as of the lLattimer & Swestvl (1199 11) (LS) 
EOS with an incompressibility at nuclear densities, K, of 180 
MeV, is employed in those simulations. On top of a strik- 
ing evidence that favors a stiffer EOS based on the nucle ar 
experimental data (K = 240 ± 20 MeV. IShlomo et al.l (120061) ). 
the soft EOS may not accoun t for the recently obser ved mas- 
sive neutron star of ~ 2M Q (iDemorest et alJl2010lPl . Using 
a stiffer EOS , the explosion energy may be even lower as 
inferred from iMarek & Jankal ( 120091) who did not obtain the 
neutrino-driven explosion for their model with K = 263 MeVQ- 
What is then missing furthermore ? The neutrino-driven 
mechanism would be assisted by other candidate mechanism s 
such as the acoustic mechanism (e.g.. iBurrows et all (120061)) 
or the magnetohydrodynamic mech a nism (e.g.. [Kotake et al.l 
(120041): ITaklwaki et al.l (12001 120091): IBurrows et all (l2007al): 
iGuilet et al.l (120101): lObergaulinger & Jankal (120111) : ?, see 
also iKotake et al.l (120061) for collective references therein). 
We may get the answer by taking into account new ingredi- 
ents, such as exotic phy s ics in t he core of the proton eutron star 
(e.g., iTakahara & Satol (fl988h : ISagert et al.l (I20091)'). viscous 
heating by the magnetorot ational instability dThompson et al.l 
l2005t iMasada et al.ll20llh . or energy dissipation via Alfven 
waves dSuzuki et al.ll2008l) . 

But before seeking alternative scenarios, it may be of pri- 
mary importance to investigate how the explosion criteria ex- 
tensively studied so far in 2D s imulations could or cou ld not 
be changed in 3D simulations. iNordhaus et al.l (120101) is the 
first to argue that the critical neutrino luminosity for produc- 
ing neutrino-driven explosions becomes smaller in 3D than 
2D. They employed the CASTRO code with an adaptive mesh 
refinement technique, by which unprecedentedly high resolu- 
tion 3D calculations were made possible. Since it is generally 
computationally expensive to solve the neut rino transport in 
3D, th ey employed a light-bulb scheme (e.g., janka & Miillerl 
dl996l) ) to trigger explosions, in which the heating and cool- 
ing by neutrinos are treated by a parametric manner. Since 
the light-bulb scheme can capture fundamental properties of 
neutrino-driven explosions (albeit on the qualitative grounds), 
it is one of the most prev ailing approximation a dopte d 
in recent 3D models (e.g. Ilwakami et ail (120081 120091) : 
IWongwathanarat et al.l (120101) ). A number of important find- 
ings have been reported recently in these simulations, includ- 
ing a potential role of non-axisymme t ric SA S I flows in gen- 
erating spins (IWongwathanarat et al.l (120101): iRantsiou et al 
(201C), see also iBlondin & Mezzacappal (120071) : iFernandez 



d2010l) ) and magnetic fields dEndeve et al.l 120101) of pulsars, 
stochas t ic nat u re of gravita t ional-w ave (e.g., IKotake et al.l 
d2009bl 1201 II): IMuller et al I d2011l) ) and neutrino emission 
(e.g.. lDuan & Knelleri (120091) ). 

To go up the ladders beyond the light-bulb scheme, we ex- 
plore in this study possible 3D effects in the supernova mecha- 
nism by performing 3D, multigroup, radiation-hydrodynamic 
core-collapse simulations. For the multigroup transport, the 
IDSA scheme is implemented, which can be done rather 



1 The maximum mass for the LSI 80 EO S is about 1.8Mq (e.g., 
IQ'Connor &TM <I(Jil) ; lKiuchi & Kotake! f2pQ8h >. 

2 On the other hand, they obtained 2D explosions for Shen EOS (K = 
281MeV, H.-T. Janka, private communication). 



in a straightforward ma nner by extending our 2D modules 
dSuwa et al.1 120101 l20llb to 3D. This can be made possible 
because we apply the so-called ray-by-ray approach (e.g., 
iBuras et ail d2006l) ) in which the neutrino transport is solved 
along a given radial direction assuming that the hydrodynamic 
medium for the direction is spherically symmetric. From a 
technical point of view, it is worth mentioning that the ray- 
by-ray treatment is highly efficient in parallizatiorflon present 
supercomputers, most of which employ the message-passing- 
interface (MPI) routines. We focus he re on the evolution of an 
1 1 .2M Q star of iWooslev et al.l d2002l) . We first choose such a 
lighter progenitor star, no t only because we follow a tradi- 
tion in 2D literatures (e.g.. lBuras et alj ([2006): B urrows et alj 
d2006l) ). but also because the neutrino-driven shock revival 
for the progenitor was reported to occur r ather earlier after 
bounce in 2D models by lBuras et al.l (12006ft . So we anticipate 
that the cost of 3D simulations would not be too expensive 
for the progenitor. By comparing with our ID and 2D results, 
we study how the increasing multi-dimensionality could af- 
fect the postbounce supernova dynamics. 

The paper opens with descriptions of the initial models and 
the numerical methods (Section 2). The main results are 
shown in Section 3. We summarize our results and discuss 
their implications in Section 4. 

2. NUMERICAL METHODS AND INITIAL MODELS 

The basic evolution equations for our 3D simulations are 
written as, 

dp 



^ + P V.v = 0, 
dv 

p— =-VP-pV$, 
df 



(1) 

(2) 



dl 



+ V- [{e*+P)y\ = - P v-V$+Q E , (3) 



= r A 



df 

A$ = 47rGp, 



(4) 
(5) 



where /9,v,P,v,e*,<J>, are density, fluid velocity, gas pressure 
including the radiation pressure of neutrinos, total energy den- 
sity, gravitational potential, respectively. 4 denotes the La- 
grangian derivati ve. As for the hydr o-solver, we employ the 
ZEUS-MP code (lHaves et al.l 120061) w hich has been modi- 
fied f or core-collapse simulations (e.g.. Ilwakami et aLll2008l 
120091) . Qe and (in Equations (0 and (HJl) represent the 
change of energy and electron fraction (Y e ) due to the inter- 
actions with neutrinos. To estimate these quanti ties, we em- 
ploy the IDSA scheme dLiebendorfer et alj|2009l) . The IDSA 
scheme splits the neutrino distribution into two components, 
both of which are solved using separate numerical techniques. 
Although the current IDSA scheme does not yet include heavy 
lepton neutrinos (y x ) and the inelastic neutrino scattering with 
electrons, these simplifications save a significant amount of 
computatio nal time compared to the canonical Boltzmann 
solvers (see [Liebendorfer et al] d2009l) for more details). As 
already mentioned, we employ the ray-by-ray approximation, 
by which the 3D radiation transport is reduced essentially to 
the ID problem. Following the prescription in IMuller et alj 
(I2010h . we improve the accuracy of the total energy conser- 
vation by using a conservation form in Equation 0, instead of 

3 along each radial ray 



FIG. 1. — Three dimensional plots of entropy per baryon (left panels) and logarithmic density (right panels, in unit of g/cm 3 ) for three snapshots (top; t = 15 
ms, middle; t = 65 ms, and bottom; t = 125 ms measured after bounce (; = 0)) of our 3D model. In the right panels, velocities are indicated by arrows. The 
contours on the cross sections in the x = (back right), y = (back bottom), and z = (back left) planes are, respectively projected on the sidewalls of the graphs. 
For each snapshot, the linear scale is indicated along the axis in unit of km. 
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FIG. 2. — Same as Figure 1 but for the net neutrino heating rate (left panels, logarithmic in unit of erg/cm 3 /s and r ac | v /7i lcat (right panels, see text for details), 
which is the ratio of the advection to the neutrino heating timescale. The gain region (colored by red in the top left panel) is shown to be formed at around t — 65 
ms after bounce, which coincides with the epoch approximately when the neutrino-driven shock revival initiates in our 3D model. The condition of T a( j v /Ti lcat > 1 
is satisfied behind the aspherical shock, the low-mode deformation of which is characterized by the SASI (bottom right panel). 
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solving the evolution of internal energy as originally designed 
in the ZEUS code. A Poisson equation (in Equation (5)) can 
be solved either by the ICCG0 method in the original ZEUS- 
MP code or by the multi-domain spectral method developed 
in the Lorene code dGrandclement & Novakll2009t) . For the 
calculations presented here, the monopole approximation is 
employed. 

The computational grid is comprised of 300 logarithmically 
spaced, radial zones to cover from the center up to 5000 km 
and 64 polar (6) and 32 azimuthal (0) uniform mesh points 
for our 3D model, which are used to cover the whole solid 
angle. To vary numerical resolutions, we run one more 3D 
model that has one-half of the mesh numbers in the <f> direc- 
tion («^=16), while fixing the mesh numbers in other direc- 
tions. Both in 2D and 3D models, we take the same mesh 
numbers in the polar direction («g=64), so that we could see 
how the dynamics could change due to the additional degree 
of freedom in the <fi direction. For the spectral transport, we 
use 20 logarithmically spaced energy bins reaching from 3 to 
300 MeV. For our non-rotating progenitor, the dynamics of 
collapsing iron core proceeds totally spherically till the stall 
of the bounce shock. To save the computational time, we start 
our 2D and 3D simulation by remapping the ID data after the 
stall of the bounce shock to the multi-D grids. To induce non- 
spherical instability, we add random velocity perturbations at 
less than 1 % of the unperturbed radial velocity. 

3. RESULTS 

In the following (section 13. U . we first outline hydro dy- 
namic features in our 3D model. Then in sections 13.21 and 
3.3, we move on to discuss how 3D effects impact on the ex- 
plosion dynamics by comparing with the ID and 2D results. 

3.1. 3D dynamics from core-collapse through postbounce 
turbulence till explosion 

Figure Q] shows three snapshots, which are helpful to char- 
acterize hydrodynamic features in the 3D model. Top panel 
is for t = 15 ms after bounce, showing that the bounce shock 
stalls (indicated by inward arrows in the top right panel) at a 
radius of 150 km. Note that colors of the velocity arrows are 
taken to change from yellow to red as the absolute values be- 
come larger. By looking carefully at the top right panel, mat- 
ter flows in supersonically (indicated by reddish arrows) in the 
standing shock (the central transparent sphere), and then ad- 
vects subsonically (indicated by yellowish arrows) to the pro- 
toneutron star (PNS or the unshocked core, the central bluish 
region in the top left panel). As seen, the entropy (left panel) 
and density (rightpanel) configurations are essentially spher- 
ical at this epochQ 

The middle panels shows an epoch (t = 65 ms) when the 
neutrino-driven convection is already active. From the right 
panel, turbulent motions can be seen (arrows in random di- 
rections) inside the standing shock, which is indicated by the 
boundary between red and yellow arrows. The entropy behind 
the standing shock becomes high by the neutrino-heating (red- 
dish regions in the left panel). The size of the neutrino-heated 
hot bubble becomes larger in a non-axisymmetric way later 
on, which is indicated by smaller structures encompassed by 
the stalled shock (i.e., inside the central greenish sphere in the 
left panel). 

4 Incomplete Cholesky Conjugate Gradient 

5 It is approximately 5 ms after we remap the ID data to the 3D grids. 



The bottom panels (f = 125 ms) show the epoch when the 
revived shock is expanding aspherically, which is indicated 
by the outgoing yellowish arrows in the right panel. The as- 
phericity of the expanding shocks could be more clearly visi- 
ble by the sidewall panels. From the entropy distribution (left 
panel), the expanding shock is shown to touch a radius of 
~ 500 km (the projected back bottom panel). Inside the ex- 
panding shock (enclosed by the greenish membrane in the left 
panel), the bumpy structures of the hot bubbles are seen. In 
contrast to these smaller asphericities, the deformation of the 
shock surface is mild, whic h is a consequence of the SASI as 
will be discussed in section l3~2l 

Figure [2] shows the net neutrino heating rate (left panels) 
and the ratio of the residency timescale to the neutrino-heating 
timescale (right panels) for the two snapshots in Figure [T] (at 
t = 65 ms (top panels) and t = 125 ms (bottom panels)). Here 
the residency timescale and the neutrino heating timescal^l 
are locally defined as, 



t KS (r,9,(f>) = 



fheatO,6>,</>) : 



r s hock(^,0)-r 



|<?bind 



for v r < 0, 
for v r > 0, 



Ql. 



(6) 



(7) 



where r ga ; n is the gain radius that depends on 9 and <fi, r s h oc k 
is the shock radius, and v r is the radial velocity. We take 
the above criteria in order to estimate the residency timescale 
for material with positive radial velocities (v r > 0) behind 
the shockQ The heating timescale can be rather straightfor- 
wardly defined by dividing the local binding energy (ebind = 
ip\ 2 + e- p$ [erg/cm 3 ]) in the gain layer by the net neutrino- 
heating rate (:g+ total [erg/cm 3 /s]). 




60 80 100 120 140 160 180 
Radius [km] 

FIG. 3. — Ratio of sound speed (c s ) and escape velocity (vesc) squared 
as a function of radius in our 3D model. This snapshot roughly coincides 
with the epoch when the neutrino-heating begins to revive the stalled bounce 
into explosion. This result is in agreement with the ante-sonic condition 
( c .i/ v Ssc ~ 0.2) for producing explosions (Peicha & ThomDso ri2011f) . 

At t = 65 ms after bounce (top panels), the gain region 
is clearly formed (reddish region in the left panel), where 

6 Originally the ratio of the advection to the neutrino-heating timescale is 
known as a useful quantity to diagnose the success {t^/t^^ > 1, i.e., the 
neutrino-heating timescale is shorter than the advection timescale of material 
in th e gain region) or failure (r a i ) v / Th ea t < 1) of the neutrino-driven exp losion 
(e.g.. lBurrows & Goshy| jl993); Janka 12001); Thompson et al. I20 0j)). 

7 W e take the word of "residency" timescale from Murphv & Burrows! 
120081) . which we feel more appropriate than the canonical "advection" 
timescale especially when we need to estimate the timescale regarding the 
postshock material with positive radial velocities. 




FIG. 4. — Top panels show distributions of the pressure perturbation (Ap/(p), see text for definition). Left and right panel corresponds to a partial cutaway 
in the FZ and XY plane, respectively. Bottom panels show the vorticity distributions ((V X v)j_ with _L being <fi or 8 in the left and right panel, respectively). 
The circle that screens between the regions colored by blue and red and the whitish region outside corresponds to the surface of the stalled shock. Note that the 
positive and negative values are colored by red and blue, respectively (e.g., +| Ap/ (p) | (red) or —\Ap/(p) | (blue) for the pressure perturbation). The linear scale 
and the time of these snapshots of the 3D model (t = 60 ms after bounce) are indicated in the top right, and bottom right edge of each plot. 



the ratio of the two timescales exceeds unity (yellowish re- 
gion in the right panel). At t = 125 ms after bounce (bottom 
panels), the ratio reaches about 2 (reddish region in the bot- 
tom right panel) behind the shock (compare the bottom right 
panel in Figure[T|i, which presents an evidence that the shock- 
revival is driven by the neut rino-heating mechanism. Recently 
iPejcha & Thompson! (1201 lb proposed an alternative definition 
of the onset time of the explosion, which is the so-called ante- 
sonic condition. From Figure [3] it can be seen that the crite- 
ria that the ratio of sound speed and escape velocity squared 
~ 0.2, is also satisfied in our 3D model when the neutrino- 
driven explosion sets in. 

Figure 21 shows distributions of pressure perturbation (top) 
and vorticity (bottom) at t = 60 ms after bounce. Here the 
pressure perturbation is estimated by Ap/(p), with Ap rep- 
resenting the deviation from the angle average pressure (: (p)) 
at a given position. Here we define the angle average of vari- 
able A as 

(A)= J 



4tt 



(8) 



The positive and negative deviations are colored by red 
and blue, respectively (e.g., +\Ap/(p)\ (red) or -\Ap/(p)\ 
(blue)). The left and right panels are for an equatorial (8 = 
7r/2, and (f> = 0) and a polar observer (9 = 0), respectively. In 



each plot, the circle that screens between the colored region 
and the whitish region outside corresponds to the surface of 
the stalled shock. From the top left panel, it is shown that the 
pressure waves (colored by red or blue) propagate outwards 
up to behind the stalled shock in a concentric fashion. Seen 
from the polar direction (top right), the color pattern at this 
snapshot indicates the dominance of £ = 2 and m = 2 modes 
in the pressure perturbation, whic h is r elated to the growth of 
SASI as we will discuss in section [3~l2l 

The vorticity distributions seen from the equator (bottom 
left panel) show that the red and blue stripes appear alterna- 
tively behind the stalled shock. Seen from the pole (bottom 
right), the vorticity waves are shown to be spinning around the 
polar axis (the origin of the figure), which would be related 
to the growth of the spiral SASI modes. These fundamen- 
tal features of the acoustic-vorticity feedbacks are akin to the 
ones obtained in ISato et af] (120091) who studied extensively 
the properties of SASI by their idealized numerical simula- 
tions. Our results might provide a supporting evidence tha t 
the advective-acoustic cycle (e.g.. lFoglizzo & Tagger! (120001) ; 
izzol j200U 120091) ) does work also in 3D simulations. 

Figure shows the spacetime diagrams of entropy disper- 
sion (a s ) for the 3D model. Note that the dispersion of quan- 
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tity A with respect to angular variation is defined by 



a A = 




(A))7(47r), 



(9) 



where (A) represents the angle average (Equation (8)). It is 
rather uncertain where the entropy production actually takes 
place in the supern ova core in the co ntext of the advective- 
acoustic cycle (e.g.. lSato et al.l (120091) V The primary position 
is the surface of the PNS, where the advecting material re- 
ceives faster deceleration by the walls of the PNS due to th e 
localized gravitational potential (e.g., iBlondin et al.l d2003l) ). 
In addition, the infalling material could also receive faster 
deceleration just outside the gain radius, where the neutrino- 
heating becomes maximum. Our 3D results tell that both of 
the two candidates are relevant indeed. As seen from Fig- 
ure 5, the position where the entropy production takes place, 
roughly coincides with the gain radius (the dotted grey line) 
before 125 ms after bounce (indicated by the upward arrow). 
Later on, the position is shown to transit to the surface of the 
PNS surface (the dotted black line). 

Until now, we have focused on the postbounce dynamics 
only for our 3D model. From the next sections, we move on 
to look into more detail how they differ from the ID and 2D 
results. 
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FIG. 5. — Entropy dispersion (cr,) in spacetime diagrams for the 3D model 
(see text for the definition). The dotted gray line represents the position of 
the gain radius, while the dotted black line shows the position of the PNS 
surface. The black arrow inserted at around 125 ms represents the epoch 
when the position where the entropy production takes place, shifts from the 
gain radius (dotted gray line) to the PNS surface (dotted black line). 



of paramet ric 3D explosion models by using the same hydro- 
code (e.g.. liwakami et all d2008ll2009T) : lKotake et ail d2009bl 
120111) ). To clearly witness the stochastic natures concerning 
the explosion direction, we may need to investigate a number 
of 3D models by changing initial perturbations and numerical 
resolutions systematically, which we think as an important ex- 
tension of this study (Takiwaki et al. in preparation). 

The left panel of Figure [7] shows mass-shell trajectories for 
the 3D (red lines) and ID model (green line), respectively. 
At around 300 ms after bounce, the average shock radius for 
the 3D model exceeds 1000 km in radius. On the other hand, 
an explosion is not obtained for t he ID model, which is in 
agreement with Bur as et al.l ( 120061) . The right panel of Figure 
[7] shows a comparison of the average shock radius vs. post- 
bounce time. In the 2D model, the shock expands rather con- 
tinuously after bounc e. This trend is qu alitatively consistent 
with the 2D result by iBuras et al.l (120061) (see their Figure 15 
for model si 12_128_f), however the average shock of our 2D 
model expands much faster than theirs. We suspect that all of 
the neglected effects in this work including general relativis- 
tic effects, inelastic neutrino-electron scattering, and cooling 
by heavy-lepton neutrinos, could give a more optimistic con- 
dition to produce explosions. Apparently these ingredients 
should be appropriately implemented, which we hope to be 
practicable in the next-generation 3D simulations. 

Comparing the shock evolution between our 2D (green line 
in the right panel of Figure [7|i and 3D model (red line), the 
shock is shown to expand much faster for 2D. The pink line 
labeled by "3D low" is for the low resolution 3D model, in 
which the mesh numbers are taken to be half of the stan- 
dard model (see Section 2). Comparing with our standard 
3D model (red line), the shock expansion becomes less en- 
ergetic for the low resolution model (later than ~ 150 ms). 
Above results indicate that explosions are easiest to obtain 
in 2D, followed in order by 3D, and 3D (low ). At first sight, 
this m ay look contradicted with the finding by lNordhaus et alJ 
( 120101) who pointed out that explosions could be more easily 
obtained in 3D than 2D. In the following section, we proceed 
to discuss what is the reason of the discrepancy more in detail. 

3.3. Comparison between 2D and 3D 

In this section, we move on to illuminate the key differences 
between our 2D and 3 D mod els. For the purpose, we high- 
ligh t the S ASI (section ["3. 3. U an d conv ective activities (sec- 
tion !3.3.2T i, the resi dency (section l3.3.3b and neutrino-heating 
timescales (section l3.3.4l i. respectively. 



3.2. Blast morphology and explosion dynamics 

Figure|6]shows the blast morphology for our 3D (left panel), 
2D (middle), and ID (right) model, respectively. In the 2D 
model (middle panel), the morphology is symmetric around 
the coordinate symmetry axis0 In contrast, non-axisymmetric 
structures are clearly shown in the 3D model (left panel). The 
direction of explosion is rather closely aligned with the polar 
axis in the 3D model. Owing to the use of the spherical co- 
ordinates, we cannot omit the possibility that the polar axis 
still gives a special direction in our 3D simulations. However 
we suspect that the alignment might be just an accident, be- 
cause the axis-free 3D explosions were obtained in a number 

8 Note that the polar axis is tilted (about it /4) both in the left and middle 
panel. 



3.3.1. SASI activities in 2D and 3D 

To compare the SASI activities in 2D and 3D, we first per- 
form the mode analysis of the shock wave. The deformation 
of the shock surface can be expanded as a linear combination 
of the spherical harmonics components Yi„,(9, (j>): 

oo / 

Rs(0^) = J2J2 c ''" Y ''" {d ^^ 

1=0 m=-I 

where Y/,„ is expressed by the associated Legendre polynomial 
Pi„, and a constant K/,,, given as 

Y lm = K lm P lm (cos6)e im *, 

1 21+1 (l-m)\ 
!m ~\] 4tt (l+my: 
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FIG. 6. — Volume rendering of entropy showing the blast morphology in our 3D (left), 2D (middle), and ID (right) model (at / = 
respectively. The linear scale is indicated in each panel. 
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FIG . 7 . — Time evolution of the 3D model, visualized by mass shell trajectories in thin gray lines (left panel). Thick red lines show the position of shock waves, 
noting that the maximum (top), average (middle), and the minimum (bottom) shock position are shown, respectively. The green line represents the shock position 
of the ID model. "1.30" and "1.40" indicates the mass in unit of M© enclosed inside the mass-shell. Right panel shows the evolution of average shock radii for 
the 2D (green line), and 3D (red line) models. The "3D low" (pink line) corresponds to the low resolution 3D model, in which the mesh numbers are taken to be 
half of the standard model (see Section 2). 



Here the expansion coefficients read, 



Cl„ 



2k 



d6 sin 0*5(0,^(0,0) 



(10) 



where the superscript * denotes complex conjugation. 

Figure [8] shows the time evolution of the expansion coeffi- 
cients (Equation ( flOb ) for the 3D (left panel) and 2D model 
(right panel), respectively. As can be seen, the amplitude 
of each mode grows exponentially until ~ 100-150 ms after 
bounce, which corresponds to the linear SASI phase. 

The top panels show that the mode of (£,m)=(2,Q) (green 
line) is dominant when the SASI enters to the saturation 
phase, which is common both in our 2D and 3D models. 
The epoch when the SASI shifts from the linear to non-linear 
phase is much delayed for 3D (t > 150 ms) than 2D (f > 80 
ms), which w a s also seen in the parametric 3D models by 
llwakami et alJ (120081) (e.g., their Figure 9). These transi- 
tion timescales are also consistent with iBuras et ai] (12006b 
who employed the same progenitor model as ours in their 2D 
simulations in which detailed neutrino transport was solved 
(see their Figure 22). It is also worth mentioning that the 
timescale seems rathe r insensitive to the em ployed progeni- 
tor. In fact, Figure 5 in lMarek & Jankal ((2009) shows the tran- 
sition timescale to be around 150 ms for a 15 M Q progenitor 
model. In the bottom panels of Figure [8] the saturation lev- 
els of the even modes of (£, \m\)= (4,0), (4,2), (4,4) in 3D are 
shown to become much larger than those in 2D (pink line), 
while the odd mode of (£, m)= (3,0) is much the same. 



Bas ed on the pioneering work by iHouck & Chevalier] 
(119921) . the linear grow th rate of the SASI in core-collapse 
case was presented by IScheck et al.l ((2008). They pointed 
out that the cycle efficiency (:Q) which represents how many 
times the average radius expands compared to the original po- 
sition per a unit oscillation frequency (:w osc ) of the SASI, is 
an important quantity to characterize the linear growth rate. 
From Figure 8, Q and ui~^ c in our simulation are approxi- 
mately estimated to be 2 and 25 ms, respectively. Note that 
these values ar e in agreemen t with the ones obtained in 2D 
simulations by IS check et alJ (120081) (e.g., their Figure 17). 
From the two quantities, the linear growth rate can be straight- 
forwardly estimated as exp(ln(Q)t w osc ), which is shown in 
the top panels of Figure 8 as black-dotted lines. As can be 
seen, the growth rates observed both in our 3D (top left panel 
in Figure 8) and 2D simulations (top right) are close to the 
linear growth rate, which seems a rather generic trend for the 
low-modes {£= 1 , 2) of the SASI. Note here that the normal- 
ized amplitude of the shock (the value in the vertical axis of 
Figure 8) shortly after bounce is actually so small as 10~ 6 to 
10~ 5 . Deduced from the linear growth rate above, the am- 
plification due to the SASI in the linear phase is at most ~ 20 
within 100 ms after bounce. On the other hand, the amplitudes 
do increase more than about 10 _1 /10~ 5 <~ 10 4 for the epoch. 
So the SASI is not the only agent for the shock deformation. 
In fact the amplitudes are observed to sharply increases from 
10" 5 to - 5 x 10" 3 within 10 ms after bounce, which is pre- 
dominantly driven by the Rayleigh-Taylor instability behind 
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with d<b/dr being the local gravitational acceleration. Ci is 
the Ledoux-criterion, which is given by 



C L = - ^ 



dp 

OP 



dP 

ds 



p,Y, 



— + 



dr 



FIG. 8. — Time evolution of the normalized amplitudes |c/ m /coo| in our 3D (left panels) and 2D (right panels) model, respectively. Lower and higher modes 
are selected in top and bottom panels. Note that the time in the figure is measured after bounce. The black dotted lines labeled by "exp-fit" in the top two panels 
indicate the linear growth rate of the SASI (see text for more details). 

the shock. To summarize, our results suggest that the rapid 
deformation after bounce is triggered by the Rayleigh-Taylor 
instability and the subsequent deformation with much milder 
growth rates is predominantly determined by the SASI. 

Concerning the saturation levels of the dominant mode of 
(£,m)=(2,Q) between 2D and 3D, it is slightly larger for 2D 
than 3D (top panels, green line). The dominance of this bipo- 
lar mode can be also seen in the blast morphology (Figure 6). 
The second-order mode of (£,m)=(l,Q) is shown to be much 
smaller for 3D than 2D (red lines in t he top panels in Figure 
|8}. This is qualitatively consistent with Nor dhaus et al.l (120101) 
who did not observe the dominance of the (£,m)=(l,0) mode 
in their 3D m odels. Note that this agreement might be just 
by chance. In llwakami et al.l (120081) . they observed the dom- 
inance of (£,m)=(l,0) mode in the saturation phase (see their 
Figure 12). From 3D results reported so far (llwakami et al.l 
l2008ll2009t[Wongwathanaratet al.ll20lol:lFernandezll2010l) . it 
seems almost certain that the low modes (£=1,2) are dominant 
in the 3D SASI-aided neutrino-driven explosions. However it 
might be rather uncertain which of the two modes (£=\ or 2) 
becomes dominant. This may reflect the fact that the explo- 
sion dynamics in 3D proceeds totally stochastically. 



3.3.2. Convective activities in 2D and 3D 

To discuss convective activities, we compute the Brunt- 
Vaisal a (B-V) frequency which is defined as (e.g jBuras et al.l 
(120061) ), 



op 

(12) 

with Y\ being the lepton fraction. It predicts instability in 
static layers if Cl > 0. The B-V frequency denotes the lin- 
ear growth rate of fluctuations, if it is positive (instability). If 
it is negative (stable), it denotes the negative of the oscillation 
frequency of stable modes. 

The left panel of Figure [9] shows the profile of the B-V fre- 
quency for our 3D model at 10 ms after bounce. The negative 
entropy gradient (bottom left panel) between the gain radius 
(<~ 100 km in radius) and the stalled shock (~ 160 km) makes 
there convectively unstable. The region in the vicinity of the 
PNS (10-20 km in radius, bottom right panel) has a nega- 
tive lepton gradient, which can make there convectively un- 
stable (the top right panel). However this region turns out to 
be convectively stable due to positive entropy gradient (com- 
pare bottom left panel). Both in our 2D and 3D models, the 
convectively unstable regions persist only behind the stalled 
shock triggered by the negative entropy gradient. 

Figure [TO] shows evolution of convective activities for the 
3D (left) and 2D (right) model, respectively. To measure the 
strength of convective activities, we define the anisotropic ve- 
locity as, 



(p[(Vr 



{v r )) 2 + V 2 g +vl 



))/</>>• 



(13) 



w B -v = sgn(C L )< 



CLd$ 
p dr 



(11) 



By this definition, higher anisotropy comes from greater de- 
viation in the radial motions (v r - (v r )) or larger non-radial 
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Radius [km] Radius 1km] 

FIG. 9. — Profile of the Brunt- Vaisala (B-V) frequency (left top panel) and entropy (left bottom pane) at 10 ms after bounce for our 3D model. The right panel 
shows the B-V frequency contributed only from the lepton gradient. Note that angle-averaged quantities are plotted here. The horizontal green line in the top 
panels is shown just for the reference of wb-v = 0. 
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FIG. 10. — Evolution of convective activities in our 3D (left) and 2D (right) model, respectively. The color-scale is for anisotropic velocity (see text for more 
details), in which higher anisotropy (colored by yellow to red) is related to active convective overturns. The dotted gray line represents the gain radius. The 
bottom panels are just zoom up of the top panels focusing on the central region. 

( v e, v <f>) motions. 

From the top left panel of Figure [10] the initial formation 
of convectively unstable regions is shown to be around 15 
ms after bounce (seen as a sudden formation of the non-zero 
Vaniso)- Subsequently, the convectively unstable regions are 



shown to advect to the center. At around 20 - 30 km in ra- 
dius, the anisotropic velocities are strongly suppressed (seen 
as a change from yellowish to bluish region at ~ 30 ms af- 
ter bounce) due to the stabilizing positive entropy gradient 
(see the left panel of Figure [9). As a result, the convective 
overturns are shown to persistently stay in the regions above 
the PNS (~ 20 - 30 km in radius) and below ~ 50 km in ra- 
dius. This is seen as a (horizontal) yellow stripe in the bottom 
left panel. Since the infalling velocities below the gain re- 
gion (dotted gray line) are so high that the convectively unsta- 
ble material cannot stay there for long. This may be the rea- 
son that the anisotropic velocity becomes relatively low (seen 
as greenish in the bottom left panel) between the gain radius 
(< 100 km) and the upper position of the yellow strip (~ 50 
km). These overall trends obtained in the 3D model are com- 
mon to 2D (right panels). In 2D, a more drastic overshooting 



of the convectively unstable material to the convectively sta- 
ble region is seen (compare the bottom panels for 20 - 40 ms 
after bounce). The area of the brought-in convectively unsta- 
ble region (equivalently the yellowish stripe) is shown to be 
larger for 3D than 2D. Such a vigorous convective overturn 
in our 3D model becomes essential in analyzing the neutrino- 
heating timescales later in section [J. 3. 41 

Having referred to the SASI and convective activities in 2D 
and 3D, we are now ready to perform analysis of the residency 
and neutrino-heating timescales. First of all, we discuss the 
residency timescale in the next section. 

3.3.3. Residency timescales in 2D and 3D 

FigureQT]depicts the streamlines of tracer particles advect- 
ing from the outer boundary of the computational domain 
through the shock wave down to the PNS. The number of the 
tracer particles that we actually injected is <~ 10 6 , however 
only the trajectories of selected particles are shown in Fig- 
ure QT| (not to make the figure filled with particles). As seen 
from the left panel, the tracer particles first go down to the 
shock wave, which is shown by the radial straight lines. Later 
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FIG. 1 1 . — Streamlines of selected tracer particles advecting through the shock wave to the PNS for the 3D model. As in Figure|4] the left and right panel is for 
the equatorial and polar observer, respectively. Each panel shows several surfaces of constant entropy marking the position of the shock wave (greenish outside) 
and the PNS (indicated by the central sphere). The linear scale is given in the right bottom edge of each panel. 
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FIG. 12. — Comparison between our 2D and 3D model at 100 ms after bounce, showing the number of tracer-particles travelling in the gain region as a function 
of their individual residency time (left panel, see text for details), and the average advection timescale as a function of the postbounce time (right panel). 

At around ~ 70 ms after bounce, the revived shock wave 
has already reached at a radius of ~ 400 km for 2D and ~ 320 
km for 3D, respectively (see the right panel of Figure |7). In 
such a shock expansion phase, the definition of the "advection 
timescale" would become rather vague. For example, the ad- 
vection timescale is longer for 2D than 3D at t = 100 ms (right 
panel of Figure [121), however this simply reflects larger shock 
radii for 2D than 3D (e.g., right panel of Figure [7)0 Also in 
the above residency-time analysis, the longer residency time 
for 2D can be seen around f res = 70-80 ms in the left panel of 
Figure [12] (seen as a dominance of the green line over the red 
line). If the onset time of explosion could be much delaye d 
after bounce (such as ~ 600 ms as in lMarek & Jankal ( 120091) ). 
the advection(or residency)-timescale analysis between 2D 
and 3D could have been made clearer in the long-lasting bub- 
bling phase. In order to see the 3D effects more clearly, we 
plan to employ a more massive progenitor (such as 15M Q ) as 
a follow-up of this study. 

3.3.4. Neutrino-heating timescales in 2D and 3D 

Now we compare the neutrino-heating timescales between 
our 2D and 3D models. From the left panel of Figure[T3j the 
heating timescale is shown to be longer for 3D (red line). As 
seen from the right panel, this is because the total net rate of 
the heating rate is generally smaller for our 3D model. To 

' In some sense artificially, this makes the advection timescale longer by 
the larger distances between R s and R g , however this does not evidently mean 
that the 2D model can gain much more efficient neutrino-heating. 



on, as indicated by the tangled streamlines, they experience 
turbulence before falling to the PNS. The low-modes (here, 
1 = 2) oscillation of the acc retion shock due to SASI activi- 
ties (as discussed in section [3.3.11 e.g., right panel in Figure 
ITTb make the residency timescales much longer for multi-D 
models than ID. If the right panel were for 2D models, the 
streamlines would be seen as a superposition of circles with 
different diameters. In contrast, the non-axisymmetric matter 
motions can be clearly seen, which is a genuine 3D feature. 

The left panel of Figure [12] compares the number of the 
tracer particles vs. their individual residency timescales be- 
tween the 2D and 3D model (for the same snapshot in Figure 
ITTb . As seen, the maximum residency time is longer for 3D 
(fres ~ 92 ms) than 2D (~ 80 ms), which is most likely to be 
the outcome of the non-axisymmetric matter motions in the 
gain region. As well known, the longer residency time is good 
for producing neutrino-driven explosions because of the long 
exposure to the neutrino heating in the gain region. 

The right panel of Figure Q~2] shows the comparison of 
the advection timescale, co nventionally emp l oyed in litera- 
tures (e.g., Equation (4) in iMarek & Jankal d2009l) . that is 
(R s -R g )/\(v r )\ with R s , R g , and (v r ) representing the angle 
average shock radius, gain radius, and postshock radial veloc- 
ity, respectively. Against our anticipation, the averaged ad- 
vection timescale is not always longer for 3D. Before ~ 70 
ms after bounce, our 3D model (red line) has generally longer 
advection timescales. However the advection timescale later 
on can be longer for 2D. 
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FIG. 14. — Same as Figure [PU but for luminosities (left) and mean energies (right) of radiated electron (v e ) or anti-electron (v e ) type neutrinos. 
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FIG. 15. — Important quantities for analyzing the properties of neutrino emission in multi-D models (see text for more details). Top left panel compares the 
velocity dispersion (top) and the average radial velocity (bottom) between 2D and 3D model. Top right panel shows the neutrino cooling rate (top) and its 
dispersion (bottom). The bottom panel compares the profiles of maximum and minimum entropy, which is indicated by ((s) +a s ) and ((.v) — tr,), respectively. 
These panels are for 50 ms after bounce. 
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understand this feature, we analyze the neutrino luminosities 
(L u ) and mean energies ((e„)), since the neutrino heating rate 
can be symbolically expressed as oc L l/ (e 2 v ) (e.g., Equation 
(23) in lJankal (120011) ). 

From the left panel of Figure [14] the neutrino luminosities 
regardless of electron or anti-electron type are shown to be 
generally larger for 3D than 2D. On the other hand, the mean 
neutrino energies are lower for 3D (right panel). Although 
the higher neutrino luminosity is advantageous for producing 
neutrino-driven explosions, the lower neutrino energies pre- 
dominantly make the heating rate smaller, thus leading to the 
longer heating timescale in our 3D model. 

The higher neutrino luminosity in 3D is du e to th e stronger 
convective activities as discussed in section [3.3.21 The left 
panel of Figure [15] compares the velocity dispersion (top 
panel) and the average radial velocity (bottom panel) between 
2D and 3D. FigureQ2]is at 50 ms after bounce, when the con- 
vection and SASI are actively in operation. 

The top left panel in Figure [15] shows that convective mo- 
tions are much more vigorous for 3D in a radius of 30-50 
km (see also the yellowish stripe in Figure IToT>. The top right 
panel shows that the neutrino cooling rate (top) as well as its 
dispersion there (erg, bottom panel) is larger for 3D than 2D. 
For 3D models computed in this work, these features are gen- 
erally maintained before the revived shock expands further 
out (typically ~ 100 ms after bounce). The bottom panel of 
Figure [15] shows that the entropy above 30 km is generally 
smaller for 3D (red line) than 2D. This is as a consequence 
of the weaker neutrino-heating in 3D than 2D. For the time 
snapshot in Figure [15] (at 50 ms after bounce), the position of 
the electron (energy-averaged) neutrinosphere is about 75 km. 
So the convection deep below the neutrinosphere (30-50 km 
in radius) is the agent to affect the neutrino luminosity and the 
mean neutrino energy. T his trend is a kin to the one observed 
in the 2D simulations bv lBuras et al.l ( 120061) . 

In Figure [16] we proceed to perform a more detailed anal- 
ysis on matter mixing behind the shock and its impact on the 
emergent neutrino luminosity. Top panels show radial veloc- 
ities for our 3D (left) and 2D (right) model within a radius 
of 100 km, in which downflows and upflows are colored by 
blue and red, respectively. The central whitish regions corre- 
spond to the PNS, which is convectively stable (hence, with 
small radial veloci ties) due to the positive entropy gradient 
(e.g., section [332] and Figure ITOt. In the vicinity of the PNS, 
downflows and upflows are visible near the equator and pole, 
respectively in the 3D model (e.g., back left and back right 
panels in Figure [16] (top left)). The bottom left panel is same 
as the top left panel but for the normalized angle variation of 
8Yy. Here we define the normalized angle variation of quan- 
tity A as 

SA=(A-(A))/a A . (14) 

We focus on anti-electron neutrino (p e ), because the luminos- 
ity of v e dominates over that of v e during the simulation time 
(e.g., left panel in Figure [l4l. Comparing the top left to the 
bottom left panel in Figure [16] it can be seen that the positive 
sign of 5Yv (reddish region in the bottom left panel) tends to 
have a correlation with the downflows (bluish region in the 
top left panel), which is vice versa for the upflows. This is 
because material with larger Yp e in the outer layers is mixed 
down to the vicinity of the PNS that possesses smaller Yp e due 
to convective overturns. This can be a possible explanation of 
the correlation between the gain (loss) in Yp e and the upflows 
(downflows) to the PNS. Note that this relation is also visible 



in our 2D model (right panels). 

Figure [17] depicts angular variations of the flow patterns in 
Figure [16] The Mollweide projection (or 47r-map) of various 
quantities is taken at a radius of 50 km. From the top left 
panel (Sv r ), downflows are shown to flow from the pole (col- 
ored by blue), while upflows are rather uniformly distributed 
in the equator (seen like a horizontal red belt). From the bot- 
tom left panel, the color pattern of blue and red reverses with 
that of the top left panel. As already mentioned, this reflects 
the correlation between downflows (upflows) and gain (loss) 
in Yp e . Reflecting the gain or loss, the neutrino cooling rate 
{&Qv e ) has a positive correlation with Yp e (compare the top 
right and the bottom left panel). The variation in the neu- 
trino luminosity that is measured at the outermost boundary 
of the computation domain (bottom right panel) has a rough 
positive correlation with the neutrino cooling rate (top right 
panel), which may agree with one's intuition. 

Finally we show the cross correlation that we take between 
the time evolution of the mass accretion rate to the PNS and 
the neutrino luminosity (Figure \TE[. As seen, a positive cor- 
relation is commonly seen during the simulation time. This 
plot may carry a message that it is important to go beyond 
the light-bulb scheme, in which the inp ut neutrino lum i nosity 
is usually kept constant with time (e.g.. llwakami et all (120081 
120091) : lNordhaus et al.l (120101) ). To take into account the feed- 
back between the mass accretion and the neutrino luminosity, 
the spectral I PSA scheme, whic h is beyond the grey transport 
scheme (e.g.. iFryer et al.l (120021) : iFryerl (12004b '). sounds quite 
efficient in the first-generation 3D simulations. 

As suggested in the right panel of Figure [7] 3D explosions 
are more easily obtained for models with finer numerical res- 
olutions^- Our results would indicate whether the advantages 
for driving explosions mentioned above could or could not 
overwhelm the disadvantages should be tested by the next 
generation 3D simulations with much more higher numerical 
resolutions. Needless to say, the 3D results (not to mention 
2D results) should depend on the sophistication of the em- 
ployed neutrino transport scheme. Regarding the gravity, we 
should first go over the monopole approximation. This may 
not be so easy task from a technical point of view, because we 
need to implement a multigrid approach to obtain a high scal- 
ability in the MPI computing. T o go beyond the N ewtonian 
gravity is also a challenging task (M iiller et al.ll2010b . Our 3D 
results are only the very first step towards a more realistic 3D 
supernova modeling. 

4. SUMMARY AND DISCUSSION 

We have presented numerical results on 3D hydrodynamic 
core-collapse simulations of an 11 .2M Q star. By comparing 
our ID and 2D results, we have studied how the increasing 
spatial multi-dimensionality affects the postbounce supernova 
dynamics. The calculations were performed with an energy- 
dependent treatment of the neutrino transport based on the 
isotropic diffusion source approximation scheme. In agree- 
ment with previous study, our ID model does not produce 
explosions for the 11.2 M Q star, while the neutrino-driven re- 
vival of the stalled bounce shock is obtained both in 2D and 
3D models. We showed that the SASI does develop in the 
3D models, however, their saturation amplitudes are gener- 

10 To say something very solid on this respect, we apparently need to study 
the effects of numerical resolutions more in a systematic manner. But the 
result we reported here is sort of best what we can do now, which it took 
more than 4 CPU months by keeping the currently available supercomputing 
facilities at our hand running. 
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FIG. 16. — Analysis of flow patterns and matter mixing for our 3D (right) and 2D (left) model, respectively. Top panels show radial velocities, in which bluish 
and reddish regions correspond to downflows and upflows that are distinguished from their local radial velocities (negative or positive). Similar to Figure 1, the 
contours on the cross sections in the x = (back right), y = (back left), and z = (back bottom) planes are, respectively projected on the sidewalls of the graphs 
to visualize 3D structures. Bottom panels show the relative angle variation of Yp e (5Yn, see text for definition). In the regions with downfl ows (bluish in the top 
panels), the sign of 8Ya tends to be positive (colored by red in the bottom panels). All the panels are at 50 ms after bounce (same as Figure fTst . 




FIG. 17. — The 47r-maps of various quantities in Figure [16] visualized by the Mollweide projection at a radius 50 km. The top left, top right, and bottom left 
panel shows the normalized angular variation of the radial velocity (Sv r ), neutrino cooling rate (SQ), and v e fraction (57^,), respectively. The bottom right panel 
is the normalized angular variation of the (anti-electron) neutrino luminosity (8La e ), which is measured at the outer boundary of the computational domain. 
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FIG. 18. — Time evolution of cross-correlation coefficient between the 
mass accretion rate (M) to the PNS and the neutrino luminosity for our 
3D model. The coefficient between two variables of A and B is defined by 

/ dn(A-{A))(fi-(a))) 

4tt 



'".■1 



-/(a A a B ). 



ally smaller than 2D. By performing a tracer-particle analy- 
sis, we showed that the maximum residency time of material 
in the gain region is shown to be longer for 3D due to non- 
axisymmetric flow motions than 2D, which is one of advan- 
tageous aspects in 3D to obtain neutrino-driven explosions. 
Our results showed that convective matter motions below the 
gain radius become much more violent in 3D than 2D, making 
the neutrino luminosity larger for 3D. Nevertheless the emit- 
ted neutrino energies are made smaller due to the enhanced 
cooling. Our results indicated whether these advantages for 
driving 3D explosions could or could not overwhelm the dis- 
advantages is sensitive to the employed numerical resolutions. 
An encouraging finding was that the shock expansion tends to 
become more energetic for models with finer resolutions. To 
draw a robust conclusion, 3D simulations with much more 
higher numerical resolutions and also with more advanced 
treatment of neutrino transport as well as of gravity is needed. 

Finally we refer to the approximations adopted in this pa- 
per. As already mentioned, the omission of heavy lepton 
neutrinos, the inelastic neutrino scattering, and the ray-by- 
ray approach should be improved. The former two, should 
act to suppress the explosion. The ray-by-ray approach may 
lead to the overestimation of the directional dependence of 
the ne utrino anisotropies (see discussions in iMarek & Jankal 
(120091) ). Although it would be highly computationally expen- 
sive, the full-angle transport will give us the correct answer 



(e.g., lOtt et all ( 120081) : iBrandt et all d2011l) ). Our numerical 
grid in the azimuthal direction is only 32 to cover 360 degrees. 
Such a low resolution could lead to a large numerical viscos- 
ity. The numerical viscosity is expected to be large especially 
in the vicinity of the standing accretion shock, which may af- 
fect the growth of the SASI. It could also affect the growth of 
the turbulence in the postshock convectively active regions, 
which is very important to determine the success or failure of 
the neutrino-driven mechanism. To clearly see these effects 
of numerical viscosity, we need to conduct a convergence test 
in w hich a numer i cal gri dding is changed in a parametric way 
(e.g. lHanke et al.l J20ll v 

A number of exciting issues are remained to be studied in 
our 3D resul t s, such as gr a vitational-wav e signa tures (e.g., 
iKotake et all d2009al l201ll) ; iMuller et al.l j201lb). neutrino 
emission and its detectability (e.g., iKistler et al.l(l20lTh \ pos- 
sibility of 3D SASI flows ge nerating pulsar kicks and spins 



( Wongwathanaratetal 



120101). The dependen c e of pr ogen- 
itors (e.g., iBuras et al.l (120061): iBurrows et all (l2007bl) ) and 
equations of state (e.g. jMarek & Jankal d2009l) ) are important 
to be clarified in 3D computations. We are going to study 
these items one by one in the near future. 

As of July 201 1, the K supercompter in Kobe city of Japan 
is ranked as the top on the "TOP 500 list of World's Super- 
compters'H From early next year, we are fortunately allowed 
to start using the facility for our 3D supernova simulations. 
Keeping our efforts to overcome the caveats mentioned above, 
we plan to improve the numerical resolutions as much as pos- 
sible in the forthcoming run, by which we hopefully gain a 
new insight into the long-veiled explosion mechanism. 
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